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Abstract: The inversion of cosh-Hilbert transform (CHT) is one of the most crucial steps 
for single-photon emission computed tomography with uniform attenuation from truncated 
projection data. Although the uniqueness of the CHT inversion had been proved [1], there is 
no exact and analytic inverse formula so far. Several approximated inversion algorithms of the 
CHT had been developed PP[2]- In this paper, we proposed a new numerical moment-based 
inversion algorithm. 

1 Introduction 

Single-photon emission computed tomography(SPECT) is a non-invasive diagnostic technology 
in nuclear medicine. It is used to image the physiological function of various organs with the 
help of radiopharmaceutical, a biochemical molecular labeled by radioisotope. Radiopharma- 
ceutical is induced into a patient, and the emitted gamma photons, attenuated by the body, 
are measured by a gamma camera rotating around it. The task of SPECT reconstruction is to 
estimate the radiopharmaceutical distribution from the measured data. 

From the mathematical point of view, the two-dimensional (2D) SPECT reconstruction prob- 
lem is to inverse the attenuated Radon transform(aRt) of p j3j. 



where p{x) and fi{x) denote the distribution of radiopharmaceutical and the attenuation coef- 
ficient of gamma photon at a; = {xi,X2)- Further, 6 = (cos 0, sin 0) and 6^ = (— sin(/), cos0) 
denote two perpendicular vectors. Although p is a function with a compact support Q in prac- 
tice, the integral ([T]) is written over (— oo, -|-oo) for convenience. 

If /i = 0, the aRt formula ([T]) is the well-known Radon transform(Rt), which is the mathemat- 
ical basis of computed tomography(CT). We can inverse Rt by the filtering back-projection(FBP) 
formula [3] from 180° projection data due to the relation that {Rp){—s,(j) -|- vr) = (_Rp)(s,0), 
where R denotes the Radon transform. If /i 7^ 0, the filtering back-projection formula for aRt 
was firstly discovered by the author of [1] (Latter, a simpler deviation was illustrated in [5]). 

In this paper, we are interested in a special case of aRt, which assumes that the attenuation 
map fi is an uniform constant in the convex domain including the support of p. With the 
assumption that /i(x) = /io(constant) inside Q and fi{x) = outside f2, the reconstruction of p 
is equivalent to solving the exponential Radon transform(eRt) of p (See [6] [3] for details) 



Since the exponential Radon transform does not have the parity property of the Radon 
transform, i.e. (i?^(,p)(— s, -|- vr) 7^ (i?^Qp)(s,0) for fiQ 7^ 0, it is a long-standing opinion that 
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{s,<p) = Rf.p{s,<j)) = I p{s9 + te^)e-^^^'^''+^'^^''^dt, 



(1) 




(2) 
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271 data is necessary for the reconstruction of p. However, accurate reconstruction from non- 
truncated projections known only over 180° was shown to be possible for any value of /Xq [ZIIH]- 

Because of practical constraints due to the imaging hardware, scanning geometry, or ionizing 
radiation exposure, the projection data are not theoretically sufficient for exact image recon- 
struction. One of the important case is the truncated data that the projection data g{(f), s) are 
only known for a limited range of s. For example, if VL is the centered disc of radius i?, the 
projection data is truncated as long as g{(f), s) is not known for all s G [— -R, B\. 

An explicit formula was suggested in pj|, in which the author reduced the reconstruction of 
2D SPECT to solving a one-dimensional integral equation called cosh-weighted Hilbert trans- 
form(CHT) by the weight-differential backprojection of eRt. This method further was devel- 
oped further in [T] , in which the author proved a similar result of the classic Radon transform. 
Recently, this result was extended to developed the reconstruction of three-dimensional (3D) 
SPECT with uniform attenuation [10]. More importantly, this method can be used to conduct 
the reconstruction of region-of-interesting from truncated projection data [HIS]. However, there 
is no explicit formula for CHT so far. Several numerical algorithms of the CHT have been pro- 
posed 

In this paper, we propose a moment-based method to solve the CHT numerically. In this 
proposed approach, the solution of CHT is represented as the result of Tricomi formula for 
Hilbert transform followed by a correction term, which is related to the moments of the under- 
lying function. The numerical simulations show that the proposed method can solve the CHT 
efficiently for a large range of the attenuation constant /xq. 

The rest of this paper is organized as follow. We review the derivation of cosh-weighted 
Hilbert transform, and present the moment-based method in section |2j Numerical experiments 
are shown in section |3} Some discussions and conclusions are given in section |4} And the 
appendix A is devoted to some numerical techniques. 

2 Cosh-weighted Hilbert Transform and the Proposed Method 

This section consists of two parts. We firstly review the derivation of cosh-weighted Hilbert 
transform(CHT) from the exponential Radon transform [T][n]. Then, the moment-based method 
is given. 

2.1 Review of Cosh-weighted Hilbert Transform 

The basic idea is similar to differential back-projection of Radon transform (yUo = 0), but with 
a weight function. Let h{x) be defined as 



where the derivative is respect to the first variable of -R^oP, i.e., {R'^^p){s,(j)) = 0), 
and X ■ 6 denotes the standard inner product of M^. By the definition of eRt ([2]), interchanging 
the partial derivative with the integration and applying the chain rule, we can obtain 




(3) 
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Here, we assume that p is continuously differentiable. Since x = {x ■ 9)9 + {x ■ 9-^)9-^, we have 
that 

iR'^^p){x-9,^) = I e'"'''9-Vp{ix-9)9 + t'9^)dt' 
Jr 

= I e^°''9-Vp{x + {t' -x-9^)9^)dt' 
Jr 

[ e^''''9-Vp{x + t9^)dt 
Jr 



Jr t d(j) 

where we use a variable replacement of t = t' — x ■ 9 for the third equality, and the last equality 
is implied by 

^p{x + t9^) = Wpix + t9^)-it^9^) 
= -t9-Vp{x + t9^). 
Inserting ^ into ([s]) and interchanging the order of intergation, we have that 

b{x) = - pix + t9^) \ldt 

Jr t 
r pfj-ot 

= - —[p{xi,X2-t)-p{xi,X2 + t)]dt 
Jr 

p{xi,X2 - t)dt + / ——p{xi,X2 + t)dt (5) 



^ Jr 't 

_ r COsh/io(x2-X-2) . ,w , .g^ 

— z/i I pyxi,X2)ux2, l^u; 

Jr '^[^1 — ^2) 

where we use variable replacements of t and x'2 = X2 + t for the first and second terms 

of (|5|, respectively. The above integrals (|5| and ^ are understood in the sense of Cauchy 
principal value. In fact, we are only interested in the reconstruction of p at {xi,X2) G f2, a 
bounded and convex region outside which p is known to be zero. Denoting by (xi, L{xi)) and 
{xi, U{xi)) the end points of the intersection of fl with the line parallel to the a;2-axis through 
the point (xi, X2) G Q, we have 

u ^ o f^'-'"'^ cosh fio{x2- x'2) , , , 

b{xi, X2) = -271 — p{xi,X2)dx2 (7) 

Jiix,) 7r(x2 - X2) 

If n is a centered disc of radius R, we have U{xi) = —L{xi) = a/-R^ — . 

Therefore, the SPECT reconstruction problem is changed into solving the so-called finite 
cosh- weighted Hilbert transform along each vertical line with (i?^Qp)(xi, 0) and {R^gp){—Xi, n) 
known. 

For simplicity, we focus on the following standard form: reconstruct a function f{t) with its 
support supp(/) C (—1, 1) from its finite cosh-weighted Hilbert transform 

/ X cosh/ii(r — t) ^, - , , , , 

J-i 7r(r-t) 
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with c^^ = J^^ f(t) cosh(/iit)(it known. We can normalize the general form of (6 ) as the standard 
form by using the following afiine transformation 

/ii = ci/io, 
/(r) = p(xi,c + Jr), 

V(^) = —b{xi,c + dT), 

where c = ^j^d d = _ Further, we have 

c^i = \ [e-'^'^«(i?;.oP)(a;i,0) + e'^^°(i?^„p)(-xi,7r)] . 
d 



If /ii = 0, by the Tricomi formula [TT], we have 



/(t)v/r^ 



Co 



-ho{s)ds 



n j_i 7r(t — s) 

Then, we can obtain / for |t| < 1 by dividing a/1 — on both sides of n9 



(9) 



2.2 Moment-based Method 

In this subsection, we focus on the inverse of the finite cosh-weighted Hilbert transform (|8 
For the sake of reference, we introduce some notations firstly. 



Tnit) 



1 



TT 



TT 



(t - r)— V(r)t/r, 



I t — s 
1 t^' 



and the m-th moment Cm = j\ f{t)t'^dt of /. Obviously, we have that Tn{—t) = (— l)"+-'^T„(t). 
By the Taylor expansion of the function cosh nit, we have 



and 



-I /"I -I °° 2k 



dr 



,2k 



(10) 



cosh{fiit)f(t)dt 



1 



C2fe- 



(11) 



4 



Therefore, h(j{t) and Cq can be represented by /lo(^) = h^^{t) — ^ j^.^'iki't), Cq = c^^ — 



k=l 



S ^fc)l'^2A;- Applying the Tricomi formula for the finite Hilbert transform|ll|, we have 



k=l 



/c — 1 \, k — 1 / 

1 °^ .,2k 1 °o ,,2fc /•! /-l /I _ _2 

1 °o ,,2fc 1 °o ,,2fc 2fc-l 1 1 /-, _ 2 

W - E tIv^- + - E E 4.-i(-l)'QT2._,_.(t), (12) 

k=l ^ k=l ^ ''1=0 



where f^^{t) = ^(c^,, - ^^j::^hf,^{T)dT) and A™ = The equation (|l2| provides 

us a relationship between the moments of / and /^^, which tells us that we can reconstruct / 
accurately if the moments of / are known. 

Define (i„ = J^^ firstly. Multiplying i = 0, 1, 2, ■ ■ ■ , oo) on both sides of 



(12), and integrating on the interval (—1, 1), we have that 



= rf2«--B2>^^^C2fc+^P^^^2L-ir2%-i)-lC2i, (13) 

where i?2i = :^ :^f=7f '^^ = ^ /_\ sin^* acia = ' for i > 0, and i?o = 1- The second 

equality is implied by the uniform convergence of the series, and the third equality holds 
because T2k{t) is an odd function for all k by the definition of Tnjt). 



2i-l 



Similarly, multiplying J^j— ^ ( i = 1,2, ■ ■ ■ , oo) on both sides of (12), and integrating on the 
interval (—1, 1), we have that 

oo 2k ^ 

C2i-1 = f^2i-l - E ^^1)! E^2i^\^2%"-V2'-1- 

Here, we use the fact that T2k-i{t) is an even function for all k. Reformulating the equations 
(13) and (14), we can derive two infinite linear systems with respect to the odd and even 
moments of /, respectively. 

QCeven d^veni (1^) 

Pcodd = dodd, (16) 
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where Ceven = (co, C2, ■ ■ ■ 

same meanings. Further, 

Qij 
P 



C2k, 



6ii + 



(ci, C3, ■ ■ ■ , C2fe+i, ■■■), and deven and dodd have the 



/^f (2^-l) 
(2j)! (2^)!! 



2fc 

_. (2(A; + 1))! 2fc+i^2fc+i-2j, 



(17) 



00 2k 
''^J + r9i-^l 2fc-1^2(fc^j)- 



^(2A;)!- 



Here denotes the Kronecker delta. Due to the uniqueness of CHT the equation systems 
(15) (16) can be solved uniquely. 

Remark 2.1 Although we can approximate f by polynomials if the moments of f are known, 
the simulations show that the reconstructed images are dominated by distortions because of the 
computed errors and the high frequency of f . 

Therefore, we can synthesize f{t)\/l — t^ by formula (12), then divide \/l — t"^ to get / finally. 
However, in fact we are only able to get a finite number of moments of / by truncating the 
equation (15) and (16). Fortunately, we only need a few low moments in computation, because 
of the rapid decrease of \cj\ with the increase of j. 

If we truncate the first 2M + 1 terms of (10) and (11), we can obtain a finite approximation 

oif{t)VT^ 



,,2k 2fc-l 

1=0 



-1)'4,_iqT2,_i. 



(19) 



Further, we have finite approximations of the operator equations (|15j) and (16). 

PCodd ~ dodd 

where dodd = (c?i, c?3, ■ ■ ■ , d2M-i) for 1 < i, j < M and deven = {do, d2,--- 
and Codd have the same meaning as deven and dodd, respectively. Further, 



(20) 
(21) 



d2M)- Further, Cg 



M 



P 



Q^ 



k=j 

u 

+ 



2k 

/^l 4 2j-lTn2j-l 

2(fc-j)' 



(2^)! 

/if (2^-1) 
(2j)! (2z)!! 



A2j-lrp2 

^2fc-l-^2 



M 

E 

k=j 



2k 
Pi 



42j rTi2i 

(2(fc + l))! 2^+1 2fc+i-2r 



(22) 
(23) 



If the integer of M is large enough, we have that the solutions to the equations (20) and (21) 
are unique by the continuity of the spectral of operators, and we can obtain the first 2M + 1 
moments of / by solving the finite linear equation systems (20) and (21). 

Based on the discussions above, we can summarize the reconstruction algorithm of SPECT 
as the following steps: 

1. performing the weight-differential backprojection (|3| to compute the CHT along each ver- 
tical line intersected with the image domain, which is used to reducing the reconstruction 
problem into solving one- dimensional integral equations(CHT); 

2. computing /^^ by Tricomi formula, di = J^^ ^^^3^ /^^(t)(jt for z = 0, 1, 2, ■ ■ ■ , 2M further, 

and the system matrices P and Q, then solving the linear systems (15) and (16) to get 
the moments q of / for i = 0, 2 ■ ■ ■ , 2M; 

3. synthesizing /(t)-\/l — by (19), then dividing a/1 — to get f{t) for \t\ < 1. 
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3 Numerical Results 




(a) (b) 
Figure 1: Left: activity distribution (a). Right: attenuation map domain (b). 

Table 1: Definition of the ehipses forming the SPECT version of the Shcpp-Logan phantom. 



Centre(cm) 1st axis (cm) 2nd axis (cm) Polar angle(°) Intensity 



(0, 0) 


0.69 


0.92 





0.5 


(0, -0.0184) 


0.6624 


0.874 





-0.2 


(0.22, 0) 


0.31 


0.11 


72 


-0.2 


(-0.22, 0) 


0.41 


0.16 


108 


-0.2 


(0, 0.35) 


0.21 


0.25 





0.1 


(0, 0.1) 


0.046 


0.046 





0.1 


(0, -0.1) 


0.046 


0.046 





0.1 


(-0.08, -0.605) 


0.046 


0.023 





0.1 


(0, -0.605) 


0.023 


0.023 





0.1 


(0.06, -0.605) 


0.203 


0.046 





0.1 



In this section, we will validate the performance of the proposed method by numerical 
experiments. For this purpose, a SPECT version of the Shepp-Logan phantom(see figure [l](a), 
[1]) was used under various conditions. Table [l] gives a description of the ellipses that form 
this phantom. The attenuation map we used is equal to a constant inside the ellipse shown in 
figure [T]^b). Although we used different values in the ellipse for different experiments to test the 
stability of the proposed methods, we only illustrate the reconstructed images in figures |2] and 
[3] for Hi = 1.5 and 3, which are sufficiently large for the medical applications of SPECT(see [1] 
for details). 

In the experiments, perfect projection data at 1000 views sampled over [0, vr] evenly with 
400 rays are created via the eRt formula (|2]). In order to test the stability and robust of the 
proposed method about noise, the data were Poisson noised by the procedure used in [T]. The 
reconstructions were discretized in a 400 x 400 grid. 

The proposed method was performed for various computer-simulated data, including perfect, 
noisy, non-truncated and truncated data. Figure |2] displays the reconstructed images from 
perfect data, while figure [3] displays the reconstruction from noised data for /xi = 1.5 and 
3. The images in top row are reconstructed from non-truncated projection data, while these 
in the bottom row are reconstructed from the truncated data such that any measurement 
corresponding to a line not passing through the rectangular box in figure [l|^a) was discarded. 

From figures [2] and |3| we can see that the presented algorithms can carry out the reconstruc- 
tion of SPECT efficiently. As our expectation, the distortions in the reconstructions raise with 
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Figure 2: Images reconstructed from the noise-free data by the proposed algorithm. 



the increase of the attenuation coefficient fii. The profiles of the reconstructions in figure |2] at 
Xi = are displayed in figure |4j Similar to [I], we have the same observation that the noise in 
the reconstructions appears strongly spatial variant, for which the reason is the weight variant 
in (|3|(see [Ij for details). 

4 Conclusions and discussions 

In this work, a numerical method for the CHT is presented. The numerical experiments show 
that the proposed method can conduct the SPECT reconstruction for practical attenuation 
constant approximately. On the other hand, we can observe that there are great distortions 
with the increasing of attenuation constant. Though numerical methods can inverse the CHT 
approximately, great efforts are needed to seek a analytic inversion formula and more robust 
numerical methods for larger attenuation constant /ii. 
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Appendix 

In this appendix, we give the analytic formulae for the computation of Tj(t) and T- . And we 
will discuss the numerical method of di, which is very important for the accurate reconstruction. 
Computation of Tn{t): Firstly, we have [H] 

To(t) = - / dr = t. (24) 

TT t-T 



Figure 3: Images reconstructed from the noisy projection data by the proposed algorithm. 



By this result, we can compute Tn(t){n > 0) by recursive relationship 



VT S — t 



TT S — t 



TT s-t nj_^ 

= sTn-lis) - Sn-1, (25) 

where 5'„ = < „ „ , and = - f^w sin^'^(a)da = ■ 

{ B2n - B2n+2 U IS CVeU ' -k J-^ ^ ! (2n)!! 

Computation of T^^{t): Firstly, we have that T„(-t) = (-l)"T(t) by the definition of T„(t). 
Therefore, we have = 0, T2nXi = for n, > 0, and T^^'^^ = = B2k+2 because 

^o(^) = ^- By using the recursive formula (25) of Tn{t), we have 

1 rl „2fc+l 

= - / -j==T24s)ds 

= - r, J sT2n~lis))ds 

= (26) 



and 



1 s^fc 

^In+l = - / ^==T2„+i(.)rfs 

TT y_i VI - s2 

1 S^fc 

:(sT2n(s) - S2n)ds 



S 



2 



'^2n^^ ~ S2nB2k- (27) 
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Figure 4: Profiles of the images reconstructed from perfect non-truncated projection data in figure [2j 

Computation of df. In order to g et the unknown function f{t) for \t\ < 1 finally, we should 
divide \/l — on both sides of (19), which implies that tiny errors near the two ends —1 and 



1 must result in large deviations from f{t). Therefore, we have to compute Cj and di precisely 
to suppress the errors on two ends. Because the computation formula for each di can be seen 
as the ^^f^—j^ weighted integral of f^^{t)t\ we use the Gauss-Chebyshev quadrature formula [I^ 
in this paper. 



n 



(28) 



i=l 



where the Gauss-Chebyshev nodes of (28) are given by qi = cos ^t'tt. 
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